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Abstract 

The swept-field experiments on magnetic molecular solids such as Fes are studied using Monte 
Carlo simulations. A kinetic equation is developed to understand the phenomenon. It is found 
that the simulations provide a quantitatively accurate account of the experiments. The kinetic 
equation provides a similarly accurate account except at very low sweep velocities, where it fails 
modestly. This failure is due to the neglect of short-range correlations between the dipolar magnetic 
fields seen by the molecular spins. Both the simulations and the kinetic equation provide a good 
understanding of the distribution of these dipolar fields. 
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I. INTRODUCTION 



As a prototype of magnetic molecular solids, consider the one generally known as Feg. 
This material consists of molecules of [Fe802(OH) 12 (tacn) 6 ]Br 8 (H 2 0)9, in which the Fe(III) 
ions of one molecule are well separated from those of a neighboring molecule. At low 
temperatures, each molecule has a spin S = 10, a corresponding all-spin magnetic moment 
of gfiBS with g ~ 2, and an Ising-like anisotropy which translates into an energy barrier 
of 22 K [H |2] . Consider next the experiment of Wernsdorfer and Sessoli [3 J , of which a 
partial and simplified description is as follows. At low temperatures (T < 100 mK) they 
first saturate the magnetization of the sample by applying an external magnetic field H z 
along the Ising or z axis, and also apply a magnetic field H x along the hard magnetic axis 
transverse to the easy axis. They then sweep H z so as to reverse the magnetization, and 
measure the rate at which the spins reverse. Astonishingly, this rate turns out to be an 
oscillatory function of H x , even though the energy barrier and the angle between the two 
energy minimizing orientations of the spin are both monotonically decreasing functions of 
H x . The interpretation of this phenomneon is that the sweeping of H z is inducing Landau- 
Zener-Stiickelberg (LZS) transitions jH [5] between the two lowest states on opposite sides 
of the energy barrier, which take place at a rate proportional to A 2 , where A is the tunnel 
splitting between these states. The rate oscillates because A(H X ) oscillates with H x , in 
accord with theory [BHE]. 

Feg is just one of ~ 10 3 magnetic molecular solids that are now known, and which have 
been the object of much study over the last fifteen years or so. Their main characteristics are 
that the spin of one molecule is large at low temperatures, and to a good first approximation, 
one may treat the spins on different molecules as non-interacting. (A comprehensive and 
authoritative review of the entire field is contained in Ref. j2]. Shorter reviews may be found 
in Refs. P-fTT].) These solids are often known as single-molecule- magnets (SMM's), because 
many phenomena may be understood, at least qualitatively, in terms of the properties of the 
total ground state spin of a single molecule in a suitable crystal field, via an effective spin 
Hamiltonian that contains anisotropy terms reflecting the overall symmetry of the molecule 
and its local environment. The oscillatory tunnel splitting mentioned above is one such 
phenomenon. Another is that the hysteresis loops are sharply stepped, where the steps 
coincide with crossings of energy levels on opposite sides of the energy barrier [T2l [TB"] . This 
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single molecule behavior has, unsurprisingly, led to suggestions and proposals for using these 
materials in devices [H] , but they are a new class of magnetic materials and worthy of study 
in their own right for the novel phenomena they display. 

Among the various experimental tools employed to study the low temperature magneti- 
zation dynamics of such solids, the swept-field or Landau-Zener-Stiickelberg (LZS) protocol 
has proven to be one of the most fruitful. When the sweep is sufficiently rapid, the accompa- 
nying change in the magnetization can be interpreted in terms of LZS transitions as already 
mentioned, and thereby provides a measurement of the tunneling amplitude between energy 
levels on opposite sides of a barrier. Tunnel splittings measured by this technique are as 
low as 1CT 8 K in temperature units. Such low splittings are beyond the reach of any other 
method. When the sweep rate is slow, on the other hand, the interpretation is not clear-cut. 
It is essential to consider the dipole-dipole interactions between different molecules, and the 
transition is influenced by the collective dynamics of all the spins. Indeed, in this case, the 
SMM designation falls short, and a more detailed analysis is called for. 

Collective, dipole-coupled dynamics of the spins are also seen in several magnetization- 
relaxation type experiments [T6H22] . Theoretical discussions and Monte-Carlo simulations 
of these experiments have been provided by [21 1231426] , and many aspects of the experiments 
are understood. The same is not true of the swept-field experiments, and we are unaware of 
any work that addresses the collective dynamics aspects. Theories that include the dipole- 
dipole interactions in a purely static [27] or mean-fieldy [T3] way cannot explain all the 
experimental behavior, especially that at low sweep fields. It is important to have a more 
complete theory since experimentalists often interpret data in terms of the original LZS 
analysis [28H3T] . In light of the fact that the spins in molecular magnets are subject to 
strong environmental influences, especially the dipolar couplings to other molecular spins, 
and that the LZS theory is based on fully quantum mechanically coherent time evolution, it 
is not at all clear that the LZS description is at all correct a priori. Indeed, as explained in 
Ref. [15J, it is a fortunate fact that it can be used even at high sweep velocities. Wernsdorfer 
has argued [29], correctly in our view, that the LZS formula should not be used without 
first verifying that one is in the fast sweep limit, and that it can in general only yield a 
lower bound on the tunnel splitting. While the use of this formula may give one qualitative 
proof of the presence of a geometrical phase in the tunneling spectrum [S], one is on shaky 
ground if one uses the extracted tunnel splittings to build detailed models of intramolecular 
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magnetic interactions. It is clear that a recent debate about these matters [28H33] is caused 
by an incomplete understanding of how molecular magnets respond to a swept magnetic 
field. 

It is the purpose of this paper to provide a more complete analysis of the swept-field 
protocol including the collective behavior of all the spins. We study the coupled system of 
molecular spins by both Monte-Carlo simulations and by developing and solving a kinetic 
equation for the joint probability distribution of the spin (up or down) of a molecule and the 
local magnetic field seen by that spin. We compare the results of these two approaches with 
each other, and with experiments. We find that the Monte Carlo simulation agrees with the 
experimental results over almost the entire range of sweep velocities, giving us confidence 
that the physical picture underlying the simulations is correct. The kinetic equation agrees 
with the simulations up to moderately low sweep velocities, but fails at still slower velocities, 
although the failure is not as severe as for previous theoretical treatments [15] . In all cases, 
we find that a model due to Kayanuma |5j of a spin strongly coupled to a bath does a better 
job of explaining the physics than the LZS model, but it is never superior to the kinetic 
equations or the Monte Carlo simulations. The main advantage of the Kayanuma model 
is that it yields a simple formula for the magnetization reversal in terms of the tunneling 
matrix element, whereas we are unable to write a similar formula for the results of out 
kinetic equation or Monte Carlo simulations. 

We note here that a kinetic equation was written down in |23J, but this equation is formal 
in that the collision or interaction term between pairs of spins is written in terms of the two- 
site distribution of spin and local magnetic field. It is thus like the first of the equations 
of the BBGKY hierarchy, which are exact, but which cannot be solved without truncating 
the hierachy. The simplest truncation is at the level of the single-site distribution itself, and 
amounts to assuming that the two-site distribution factorizes. However, rather than obtain 
the kinetic equation in this formal way, we derive the equation by considering the various 
interaction processes which cause this distribution to change. This approach gives greater 
physical insight into the approximations made. The single-site approximation breaks down 
at ultra slow sweeps as we explain, leading to the divergence from the simulations mentioned 
above. Including two-site correlations would lead to an immensely more complex numerical 
problem, however, and experience with other problems suggests that if two-site correlations 
are indeed important, then we are faced with a humdinger, as three-site, four-site, and 
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higher-order multi-site correlations are also likely to be important. 

The rest of the paper is organized as follows. In Sec. [TT| we provide background infor- 
mation on Fe§, the LZS and Kayanuma models, and on a theoretical model for SMM's that 
incorporates the influence of the environment [15l We also describe how this model 



differs from our earlier rate equation approach [26] . In Sec. Ill , we discuss our Monte Carlo 



protocol, and the results of our simulations. The kinetic equation is discussed in Sec. IV We 
first derive this equation, and then discuss how various partial sums of the collision term in 
this equation may be physically interpreted. We also discuss how we integrate this equation 
numerically, and the results of this integration. We summarize our conclusions briefly in 
Sec. ED 



II. BACKGROUND INFORMATION AND MODELS 
A. Independent and Coherent Spin Approximation 

Let us first assume (counterfactually) that the interactions between different molecular 
spins may be ignored. The dynamics of the spin of a single molecule would then be governed 
by an anisotropy Hamiltonian of the form 

U = -k 2 S 2 z +U ± -gfi B S-H, (2.1) 

where S = (S x , S y , S z ) is the total spin of the molecule, the first term in H is the leading 
anisotropy, g is a g- factor, /i# is the Bohr magneton, H is the external magnetic field, and 
T-L± is a term in the transverse spin components that is off-diagonal in the S z basis and 
reflects the intrinsic higher-order anisotropies of the molecule. For this reason it is time- 
reversal invariant, and so contains only even powers of the spin components. In Fes, for 
example, 

Hx = (h - k 2 )S 2 x - C{S\ + Si). (2.2) 

The various parameters in T-i are well known for the most highly studied molecules. For 
Fe 8 , S = 10, g ~ 2, k x ~ 0.33 K, k 2 ~ 0.22 K, and C ~ 29 /iK. In Mn i2 , by contrast, the 
anisotropy is tetragonal based on the symmetry of the molecule, but there are believed to be 



biaxial terms of the same type as in Eq. (2.2) arising from chemical disorder, variable waters 



of crystallization, variant chemical species, etc., and there are also additional longitudinal 
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terms such as BS Z . These details are largely immaterial for this paper, since we are interested 
in ultra-low temperatures only. We may therefore focus on the lowest two energy levels, 
m = ±S, and presuppose that there is an amplitude per unit time, —iA/2H, to tunnel 
between these levels. We also assume that the effects of the transverse fields H x and H y 
have been incorporated in A, and only the longitudinal field, H z , is not. It is this field that 
is swept. 

Under these conditions, each molecular spin may be described as a two-level system 
governed by an effective single-spin Hamiltonian, 

1 (<t) A \ 

n eS = - , (2.3) 

2 ^ A -e(t) J 

where e(t) is the energy of the m = S state relative to that of the m = —S state, given by 

e(*) = 2Sgfi B H z (t). (2.4) 



We shall refer to e as the bias on the spin [35]. In writing Eq. (2.4), we have defined the 
zero of H z so that the levels ±S are degenerate at H z — 0. In the laboratory, a nonzero 
offset field may be required to cancel demagnetizing fields and bring this degeneracy about; 
H z is supposed measured from this offset. We will interchangeably refer to the two states 
either as m = ±S states, or as pseudospin-1/2 states |t) = \m — +S) and ||) = \m = —S) 
states, or as "up" and "down" states. 

We now suppose that the longitudinal field is swept at a steady rate so that 

e(t) = it, (2.5) 

and the spin is in the lower energy state, | f) as it -> — oo. Then, the probability that spin 
will flip into the | J, ) state as t — > oo is given by the classic LZS formula 

P LZS = l-exp(-7rA 2 /2|e|). (2.6) 

The limits of fast sweep, e ^> A 2 , and slow sweep, e< A 2 , are worth noting: 

e > A 2 , 



Plzs « { m (2.7) 



1, e<A 2 . 
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B. The Kayanuma model 



As mentioned in Sec. [T| the LZS model is inadequate to describe the actual experiments. 
An alternative single-spin model is that of Kayanuma [3] wherein the bias field e(t) has added 
to it a fluctuating part rj(t), which is taken as a Gaussian random process. In the limit where 
this process has very large amplitude and is delta-function correlated, corresponding to a 
very rapidly fluctuating bias, the probability that the spin will end up in the state | I) 
having started in | ^ ) is found to be 

P K = ^(l-exp(-7rA 2 /|6|)). (2.8) 

The fast and low sweep limits of the reversal probability are now given by 

7TA 2 

P K * <! 2 ^ VIS)) 




This is identical to the LZS formula for fast sweeps, but very different for slow ones. Whereas 
the LZS process describes coherent adiabatic reversal of the spin, the Kayanuma process 
describes a spin which is able to make many transitions between the up and down states 
and is therefore completely randomized at t = oo. We may think of the fluctuating field r)(t) 
as a qualitative way of describing the dipolar field of other spins, which are also undergoing 
transitions between up and down states. 



It was shown in Ref. JT5] that we also obtain Kayanuma's answer (2.8) if we assume that 
the external field is swept uniformly, and each spin flips between the up and down states 
independently of the others but with a bias-dependent rate. That is, if the probability for a 
spin to be in the up state is denoted p^(t), we take 

p r = T[e{t)]{l-2p r ). (2.10) 



Equation (2.8) follows if 



r(e)rfe=^A 2 , (2.11) 

and e(t) = it. Since this requirement on T(e) entails only the tunneling amplitude A, the 
details of the physical decoherence mechanism that justifies writing down a rate equation 
in the first place are irrelevant. In Ref. [IS], these mechanisms involved the environments 
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of nuclear spins and other molecular spins. The bias-dependent flip rate that we employ in 



this paper [see Eq. (2.13)] is a special case of that in [T5], and satisfies Eq. (2.11). 

The Kayanuma model is attractive because it is simple and economical. It is a much 
better model for the reversal than the LZS process. It is plausible that for a large collection 
of molecular spins where each one may flip at different times in the sweep cycle because 
the offset field is not the same for all spins, the field seen by any one spin has a stochastic 
character. However, it is clearly an approximation to assume that the time-dependence of 
the bias on any one spin is uncorrelated with the configuration of the other spins. As shown 
in [T5], it is not even enough to take account of the spin configuration in a mean-fieldy 
way through the spatially averaged demagnetization field. Therefore, we must consider the 
interaction between the spins explicitly. 

C. Environmentally influenced interacting spins 

In reality, the spins are neither isolated nor noninteracting as noted before. As discussed 
in [T5l |3l], there are two types of interactions to consider. First, each molecular spin is 
coupled to nuclear spins in its vicinity. These spins have the effect that the tunneling 
between the m = ±S states becomes fully incoherent for typical molecular magnets. The 
IT) ^ 14) probability for the zth spin is given by 

Pflip.i = dt, (2.12) 
where 

Here, W ~ 10Ed n with Ed n being the dipole-dipole interaction energy between the molecular 
spin and the nearby nuclear spins, and £j is the bias on site i (35]. For Fes, Ed n ~ 1 mK, 
and e.j ~ 0.1 K in temperature units (see below). 

Second, each molecular spin is coupled to all the other molecular spins via the dipole 
interaction. At first sight, this leads to a fully many-body quantum mechanical problem. 
Since the nuclear spins already render each molecular spin slow and incoherent, however, the 
field of one molecular spin on another may be taken as a c-number and quantum mechanical 
back action may be ignored. Let us define an Ising spin variable <7, on every site i, such that 
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Oj = ±1 corresponds to m.; = ±S. The dipolar part of q is then given by 

ei,dip = ^ KjjVj, (2.14) 

1-3^V (2.15) 

Here, _Ed m is the interaction energy scale between neighboring molecular spins, and a is 
their separation. Further, is the distance between spins i and j, and z^- is the difference 
between their z-axis coordinates. We estimate the dipolar field in Fes as ~ 100 Oe, which 
implies the energy scale ~ E^m ~ 0.1 K quoted above. 

The net result of these two couplings is as follows. The ratio e^jW is large for most 
molecules most of the time, so these spins are essentially frozen. A spin is unfrozen only 
when the field it sees is essentially zero (more precisely of order W or less) . We refer to this 
range of bias values as the reversibility region. In magnetization relaxation experiments, 
the combination of a static and narrow reversibility region leads to very slow and non- 
exponential time decay, which many workers have investigated theoretically [23H26] . If the 
external field is swept, however, the odds of a spin being able to flip are greatly increased. 
Particularly if the sweep rate is low, the complete dynamics can be very rich and complex, 
as we now discuss. 

Let us focus on one spin as the bias is being increased. Since the dipole energy is so 
large compared to all other energies in the problem, the flip of even a far away spin can 
change the bias on the first or central spin from a big negative value (relative to W) to a big 
positive one. The central spin is then shunted past the e = region and does not flip. This 
gives a qualitative explanation of why the net magnetization reversal can be less than that 
in Kayanuma's model. However, the central spin can be returned to the region of negative 
bias if another spin not too far away were to flip subsequently. Thus, the actual amount of 
spin reversal is hard to calculate, especially analytically. We therefore turn to Monte Carlo 
simulations of the LZS protocol. 



- ^ ( 

in \ 



ILL MONTE CARLO SIMULATIONS 



The physical problem we have described in Sec. |II C| is characterized by four parameters: 
the molecular spin dipole-dipole energy Ed m , the reversible region width W, the external 
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bias sweep rate e a , and the tunnel splitting A. For purposes of analysis and Monte Carlo 
simulation, however, these parameters can be reduced to just two dimensionless variables: 
the scaled reversible region width w and the scaled sweep rate v, given as 

W 

w = " (3.1) 



(3.2) 



All bias energies are measured in units of E& m . We use the symbol e (note the slightly 
different font) for the dimensionless biases, both applied and dipolar. It is not always 
optimal to frame the discussion in terms of dimensionless variables, and we shall back and 
forth between dimensionless and dimensionful descriptions as needed. 



A. Simulation protocol 

Our Monte Carlo simulation models a finite system of N spins, labeled <7j, that can be in 
either the up or the down state (<7j = ±1). The bias at site i is the sum of the externally 
applied bias field e a (t) and the dipole field e^dip: 

ei(t) =e a (t) + e lAip (t). (3.3) 

We have explicitly shown that the dipole field is time dependent since the spin state of the 
system is time dependent. The spins are taken to lie on a cubic lattice with lattice constant a, 
and the sample is taken to be spherical so that the initial bias distribution is approximately 
a delta function centered at the value e a (0). For the majority of the simulations, the number 
of spins, N, is chosen to be 7,153, corresponding to 25 spins on the sphere's diameter, and 
e a (t) is taken to vary linearly in time, 

e a (t) = e a t. (3.4) 

Though the number N is significantly less than that used for previous simulations of relax- 
ation we found that the results converged even for N > 1,419 (diameter > 15 spins). 
In simulations of relaxation, large values of N were necessary to ensure that the reversible 
region was never depleted entirely. In simulations of LZS sweeps, the reversible region moves 
across the bias distribution, so depletion is a much lesser concern. 
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In each Monte Carlo timestep (which we arbitrarily denote by dt), the spin flip proba- 
bilities are calculated for every spin in the system from the rate function T(e). To save on 



computation time, for most of our simulations, we used a simplification of Eq. (2.13) that 
is only nonzero for |e| < W: 

T'(e) = 7 ^Q(W-\e\). (3.5) 

(The constant multiplying Q(W — |e|) is chosen so that the integrals of T(e) and T'(e) over 
all e are identical, and Eq. ( |2.11[ ) holds with V.) In practice, this rate function gives very 



nearly the same results as Eq. (2.13) and allows a simpler view of the spin flip process: the 
only spins that are able to flip are those strictly within the reversible region, a window of 
width 2W around zero. After the spins have flipped, the change in bias at every lattice site 
is calculated and the simulation moves to the next timestep. 

Over the course of one run of the simulation, the external bias is swept from —25Ed m to 
25E dm . If de a is the change in the bias in one time step, the probability for a spin in the 
reversible region to flip in that timestep is 

ttA 2 , 

p ^ = ~w dt 

n E dm A 2 de a 
4 W E dm \e a \ 

= de a . (3.6) 

4wv 

Naturally, care must be taken to ensure that pm P <C 1. In practice, we required that it did 
not exceed 0.1. We also want to ensure that we do not skip past the reversible region in 
just one time step, so de a is also chosen not to exceed 0.1u>. This becomes computationally 
intensive for slow sweeps, and for v = 0.1, for example, we are only able to study values of 
w exceeding 0.01. 

The quantities that we record during the simulation are the magnetization per site, 

M = lj>, (3.7) 

i 

and the bias field distributions, f±(e), defined so that f±(e)de is the fraction of sites with 
spin a = +1 or —1 and a bias field in the range e to e + de. We are especially interested in 
the final magnetization, Mf. 

We have already touched on the feasibility of performing the simulations with spheres 
of as few as 15 spins on the diameter. However, to test the stability of the answers with 
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system size, we have performed some simulations with as many as 55 spins on the diameter. 
Further, in some simulations, we have employed a triangular wave for e(t) as in some of the 
experiments [21 [201 EE]- The excursion of the bias is again ±25E dm . In these cases, e refers 
to the value of \de/dt\ on any leg of the wave. Finally, a few simulations are done using the 



full Gaussian rate function (2.13). We shall mention these exceptional cases as they merit. 



B. Simulation Results 
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FIG. 1: (Color online) Monte Carlo results for the magnetization of the system as a function of 
the externally applied bias as it is swept from negative to positive values for various values of v 
and w = 0.05. Since e a varies linearly with time, this is essentially a plot of M vs. t. Each curve 



is an average of 100 runs for a system of N = 7,153 spins using the modified transition rate (3.5). 



The triangles on the right edge of the plot show AT/k> the final magnetization for each value of v 
as per Kayanuma's model. 

In Fig. [I] we show a plot of the time dependence of the magnetization for various values of 
v and w = 0.5. Readers will note that irrespective of v, the major change in M starts when 
e ext reaches a value close to zero, and there is little change once e ext exceeds ~ 5. Readers 
will also note that there appears to be a small drop in M before the main one. The reason 
is that in the initial state, when all the spins are up, the bias field is zero at almost every 
site. However, there are always some spins with nonzero bias on the surface of the sample, 
though their fraction becomes smaller as N increases. Thus the initial bias distribution is 
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not quite a delta function centered at zero bias, and the surface spins start to flip before the 
external bias reaches —w, and this is responsible for the small dip. 




FIG. 2: (Color online) Final Monte Carlo (MC) magnetization Mf as a function of the scaled 
sweep rate v for various values of w. Each data point is an average of 100 runs for a system of 



N = 7,153 spins using the modified transition rate (3.5). Also shown is Eq. (3.8) for M^k, the 
final magnetization in the Kayanuma model. 

As can also be seen in Fig. [TJ the final magnetization, Mf, agrees with the Kayanuma 
result only for large enough v (the agreement with the LZS result is even poorer). To make 
this point clearer, in Fig. [2] we show Mf vs. the sweep rate v for three different choices of w, 
along with the Kayanuma result 

M fiK = 1-2P K = e - wA2 /l**l = e-*/ v . (3.8) 

There is good agreement with the Kayanuma result for fast sweep rates, with v > 10. In 
fact for these velocities, the Kayanuma and LZS models both give good answers, but we 
do not compare the Monte Carlo answers with LZS because the disagreement for v < 10 is 
much worse. Below v of about 10, however, even the agreement with Kayanuma's model 
starts to degrade, and we obtain a residual magnetization even for very slow sweep rates. 
The agreement is worse for smaller values of w, with a lower value of w corresponding to a 
larger Mf. This is qualitatively understandable as larger w allows for greater opportunities 
for relaxation. What is interesting is that for v > 10, there is virtually no w dependence 
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in Mf, and the divergence in Mf starts to appear at about the same value of the scaled 
velocity, v ~ 10, as that where the disagreement with LZS and Kayanuma starts to show 
up. 
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FIG. 3: (Color online) The value of A that we would infer for Fes if we applied the LZS formula for 
Mf to the Kayanuma model, our Monte Carlo (MC) simulations, and the experimental data from 
Wernsdorfer et al. [3J (open diamonds). Since the actual dependence of Mf on v is incorrectly given 
by LZS, this leads to a sweep-rate dependent answer for the inferred A. We deduce the conversion 
factor from e to dH/dt by exploiting the fact that for high sweep rates, LZS and Kayanuma both 
give the correct answer for A, viz. 1.12 x 10 _7 K. 

We now show what value of A we would infer if, as experimenters often do, we use the 
LZS formula for the final magnetization [SB]. In other words, denoting this inferred value 
by A infj Lzs; we use the formula 

.9 2 lei /1 + M f \ , . 

Ainf,Lzs = -^ln(^-^J. (3.9) 

The results are shown in Fig. [3j We show what we get if we use the Kayanuma formula 
for Mf, the experimental data of Ref. [3J, and Mf as given by our Monte Carlo simulations 
with w = 0.01 and w = 0.1. From this plot, one cannot tell which value of w is better, 
although as Fig. [2] shows, w = 0.01 gives a better fit. Since we estimated W ~ 10 mK and 
Edm — 100 mK, w = 0.1 certainly seems reasonable, and even w = 0.01 is not out of the 
question. This graph shows that the magnetization reversal is not a very sensitive function 
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of w, and so is of limited value in trying to learn about nuclear spin decoherence. As can 
be seen, the disagreement between even the Kayanuma model and the experimental data 
is substantial (observe the logarithmic scale for v), while that between the data and our 
simulations is small. 




FIG. 4: (Color online) Bias distributions for a sample of N = 20,479 spins with v = 0.1, w = 0.05 
at (a) e a = -1.0E dm , (b) e Q = 0.0E dm , (c) e a = l-0E dm , (d) e a = 25.0E dm , all averaged over 100 
runs. The reversible region of the bias distribution is very narrow on the scale of this figure, and 
shows up as the vertical line. 

We now attempt to understand the nonzero residual final magnetization. To this end, 
we consider the bias distributions /±(e) over the course of a field sweep. Initially, the bias 
distribution at up spin sites is sharply peaked at zero, with only spins on the surface of the 
sample having a significant bias. As the external bias e a approaches zero, spins begin to 
flip, causing the distribution to widen. Even though the dipolar interaction is long-ranged, 
when a spin flips from up to down, the resulting change in bias on the majority of the 
sample is small on the scale of W. However, spins within a couple of lattice spacings of the 
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spin that flipped will experience a dramatic change in bias. Of the six nearest neighbors 
of a flipping spin, the four in the xy plane will experience a change in bias of — 4E , £ j m , and 
the two along the x axis will experience a change 8Ed m . This can be seen in Fig. |4Jd as 
the bias distribution has two minor peaks spaced approximately —4Ed m and 8Ed m from the 
central one. Suppose that we are at a point in the sweep where a particular test spin in 
the sample is seeing a net bias close to zero, and one of its nearby spins (not necessarily a 
nearest neighbor) flips. We refer to the second spin as the triggering spin. Let the first spin 
undergo a large negative change in bias as a result of this spin flip event. It will then be 
displaced in bias far from the reversible region, yet since the external bias is being swept 
from negative to positive values, it will reenter the reversible region at a later time provided 
that the dipolar contribution to the bias seen by it has not changed significantly. On the 
other hand, if the test spin undergoes a large positive change in bias, it will be displaced in 
bias far past the the reversible region. Unless the triggering spin flips yet again or another 
nearby spin flips so as to change the bias on the test spin by a large negative amount, it 
will continue to see a large positive bias, and remains in its original orientation, up. It will 
essentially have been shunted around the region of reversibility. These shunted spins give 
a net positive contribution to Mf, while unshunted spins will contribute an amount that is 
close to the Kayanuma result on average. 
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FIG. 5: (Color online) Final Monte Carlo magnetization as a function of the reversible region 



width w for several small values of v, using the Gaussian rate function (2.13). Each data point is 
an average of 100 runs for a system of N=7,153 spins. 

We conclude this subsection by discussing two atypical simulations. First, to confirm 
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that the dependence of Mf on w for small values of v was not an artifact of using the 



modified spin flip rate (3.5), we ran the simulation for fixed v < 10 and variable w using the 



original Gaussian rate function (2.13). The results are plotted in Fig. [5jand show that the 
w dependence is weak, although it is more pronounced for smaller v and Mf tends to the 
Kayanuma result with increasing w. 



10 



10"' 



10" 



0.0 0.5 



i • 


— Kayanuma 
♦ ♦ MC v=0.1 


♦ 


♦ 



1.0 1.5 2.0 

N 



10" 



10 



10"' 



10" 



\ * 


— Kayanuma 
♦ ♦ MCv = 1.0 


\ * 

\ ♦ 


♦ 

♦ 

♦ 



0.0 0.5 1.0 1.5 2.0 2.5 3.0 

JV 



10" 



10" 



10' 



10"' 





— Kayanuma - 
♦ ♦ MC u=10.0 ■ 




s. ♦ 

♦ 



10" 



10" 



10"' 



10" : 





— Kayanuma 
♦ ♦ MCv = 100.0 







V 



20 40 60 80 100 

N 



FIG. 6: (Color online) Magnetization M plotted against number of cycles ./V for an applied trian- 



gular wave using the original Gaussian flip rate (Eq. (2.13)) with w = 0.05 and (a) v = 0.1, (b) 



v = 1, (c) v = 10, (d) v = 100. The grey line shows the prediction of the Kayanuma model. 



Second, we show the results of simulations in which e(t) is a triangular wave. In Fig. [6 
we show the dependence of M on the number of cycles of the wave for different values of v. 
Note that we show M every half-cycle, where a cycle is a complete period of the wave. The 
relevant point here is that M decreases essentially exponentially with the number of cycles, 
except for very small v, showing that it is the fraction AM /Mi that is the same each time 
the reversible region is traversed, where Mj is the initial fraction before the traversal. This 
is true even when the agreement with the Kayanuma model is poor. 
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IV. KINETIC EQUATION 



In Ref. [26J, the results of Monte Carlo simulations of magnetization relaxation were 
analyzed in terms of coupled rate equations for three quantities: the magnetization M, the 
magnetization of the spins in the reversible bias region M r , and the number of spins in the 
reversible region N r . However, the equations did not form a closed system, and also involved 
a functional J 7 , given by 



The distribution /o-(e) was then handled via an approach called the Three Gaussian Ap- 
proximation (TGA), wherein it was modeled as a sum of three Gaussians, centered at biases 
of 0, corresponding to a spin with no neighbors flipped, and —4Ed m and SEd m , corresponding 
to spins with one nearest neighbor flipped. The widths of the Gaussians were assumed to 
be equal, and this common width and the heights of the three Gaussian were determined by 
matching the first three moments of Y2 a f<r( e ) wr th a model in which every spin was up or 
down with probabilities (l±M)/2 independently of the others. The distribution determined 
in this way depends only on the system's magnetization, and leads to a closed system of 
rate equations. These equations were then found to describe the short to moderate time 
behavior of the magnetization well. 

It is plain that this approach is inadequate for analyzing the swept-field problem. It is 
paramount to have a good approximation for the bias distribution near the reversible region 
to properly capture the probability of spin flips. For the relaxation problem, the reversible 
region is static and centered at zero bias, whereas in the case of swept field, the reversible 
region moves over the full range of the bias distribution. As a result, a larger number of spins 
are capable of flipping, and the bias distribution is altered over essentially it entire range. 
Thus, the TGA is fundamentally invalid and one cannot really identify just three peaks. The 
approximation may perhaps work for very fast sweeps since the fraction of spins that flip 
is then small, but even in this limited success it does not provide a physically satisfactory 
explanation. 

It is therefore necessary to analyze the time evolution of the full bias distribution / CT (e). 
This is naturally done in terms of a kinetic equation. We develop this equation below in 
somewhat more generality than is necessary for the precise problem of this paper with a 




(4.1) 



18 



view to applying it to other settings in the future. The basic process responsible for changes 
in the distributions are exactly the same as in Ref. [26J, so our present discussion is more 
brief. 



A. Derivation of kinetic equation 

Let us denote the rate at which a spin at bias e flips from o to —a by 

r^(e), (4.2) 

where we write a for —a in the suffixes. In general the rates T a ^(e) and r S(T (e) need not be 
equal, although in this paper they are. One way in which / CT (e) can change is by what might 
be called direct or one-spin processes, wherein in a time dt, a spin flips at any given site i, 
but does not undergo any change in the bias seen by it. This gives a contribution 



dt 



= -V^(e)f a (e)+V^(e)U(e). (4.3) 

direct 



Next, we consider two-spin processes, where we focus on a spin at a particular site, i, and 
allow the bias seen by this spin to change by virtue of a spin- flip at another site j . We refer 
to the second spin as the triggering spin, and to the first as the central spin. The latter 
is taken to not flip, as the probability for a process where it also flips is of order (dt) 2 and 
therefore negligible. There are then two processes which cause / CT (e) to change: 

site i site j site i site j 

I. a,e a',e' -> a, e" -a',e' (4.4) 
II. cr,e" a',e' -> a,e -a',e' 



Processes I and II can be thought of as loss and gain processes, respectively, since they lead 
to spins being knocked out of or into the bias region (e, e + de). Let us consider process I 
first. Since the spin at site j flips from cr' to —cr', the change in this spin is —2a'. In order 
for the bias at site % to change as shown, we must have 

e" = e - 2K ij( j\ (4.5) 



K l3 = \{e - e")a' . (4.6) 
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If we let the final bias on site % lie in the range (e", e" + de"), the number of triggering spins, 
i.e., the number of spins that satisfy this condition is 

<7(i(e-e>')x^, (4.7) 

where 

g{K) = ^ 5^ - K^) (4.8) 
is the density of couplings, by which we mean that g(K)dK is the number of sites which 



couple to the central site with couplings between K and K + dK. The sum in Eq. (4.8) is 
over an infinite lattice, and therefore independent of site i. The probability that a spin on 
a triggering site will indeed flip in time dt is T a i a i(e') dt. To find the net loss in / CT (e), we 
must multiply the number of triggering sites with the probability of a flip at those sites and 
the fraction of sites i and j that have the stipulated spins and biases. We must then sum 
over all possible values of a', e', and e". In this way we get 



dfa_ 
dt 



-T,j de ' J d ^^)9{W-e'')a')fAe)Ue)- (4-9) 

cr' 

The calculation for process II proceeds in identical fashion. This time we get 

n = E/ rfe 7 d ^^)9{W'-e)a')Ue)Uey (4.10) 



dfg 

dt 



Together, the last two equations describe what might be called the collision integral. 

Lastly, let us incorporate a swept or explicitly time dependent applied field. If there were 
no spin flip processes at all, then a site which had a bias e at time t, would at time t + dt 
have a bias 

e(t + dt) = e(t) + de a , (4.11) 
where de a is the change in the applied field in the interval dt, and we would have 

f a (e(t),t) = f a (e(t) + de a ,t + dt) 

= f^^ + ^dea+^dt. (4.12) 



Therefore, 

dfg df a de a 

dt sweep " de dt ' 



(4.13) 
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Adding Eqns. Ob, Ob, flOOl), and (|4~13j), we get 



dfa 
dt 



dfa de a 
de dt 



E 



-T^{e)f a (e)+T a9 (e)f 9 (e) 



de' de" 



f a ,(e')T aW (e') - e>')/ CT (e) - <?(§(e" - e y)/„(0 



(4.14) 



This is the kinetic equation we are seeking. We can rewrite this equation in various ways 
which help understand its structure, and also suggest algorithmic simplifications for numer- 
ical integration. So, we note that in the collision term, the integral over e' can be factored 
out of the expression completely, reflecting the fact that the precise value of the bias at the 
triggering site is irrelevant to whether a flip on this site alters the bias on the central site 
by a given amount; what matters for that is where the triggering site is located relative to 
the central site, and the spin on the triggering site. We may therefore define a net rate per 
site for spin flip, averaged over the bias distribution, and over the entire sample, 



T 9o {t) = J def a (e,t)T a(7 (e). 



(4.15) 



As indicated, this rate changes with time because the bias distribution changes. It is in fact 
a functional of the distribution. In terms of this rate, we may write the kinetic equation as 
dfa dfa de a 



de dt 



Taa(e)fa(e)+Taa(e)fa(e) 



J^Ta'a' / de' \ g {l{e-e')a')f (T {e)-g(\{e'-e)a')Ue') 



dt 



We have also changed the dummy variable of integration from e" to e' . 

An elementary but useful check on the kinetic equation is that the quantity 



(4.16) 



E / de f°& 



(4.17) 



should be time-independent, since it is just unity by normalization. Summing Eq. (4.16) 
over u and integrating over all e, we get 



d 
dt 



E / de f^ 



dt 



dfa 

de 



/ de [Taa(e)fa(e) - V c - C {e)h{e) 

a 

J^Ta'a* 1 1 dede' \g(l(e - e')a')fa(e) - g(\{e' - e)a')fa(e') 

(4.18) 
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The sweep term integrates to zero directly, since f a {e) must vanish for e — > ±00. The direct 
terms also add to zero. To see that, we first note that a sum over a is equivalent to a 
sum over a. Recognizing this, and writing the second direct term as a sum over a, and 
then interchanging the index labels a and a, the two direct terms are seen to cancel each 
other identically. To see that the two-spin terms also add to zero, we simply interchange the 
integration variables e and e' in the very last term. The two parts of the collision integral 
are then identical and they cancel each other. Hence, 



d 
dt 



£ / deU{e) 



0, (4.19) 



as desired. 

At this point it is useful to discuss the nature of g(K), the density of dipole couplings. 
As shown in Ref. [To] . 

That is, for small K, g(K) is approximately an even function, and diverges as 1/K 2 . The 



integrand of the e' integral in Eq. (4.16) therefore has a singularity at e' = e of the form 



(e' - ef 

which is equivalent to 



(4.21) 



1 dU{e) (4.22) 



e' — e de 

This is an integrable singularity if we view the integral over e' as a principal value. Hence, 
our kinetic equation is mathematically well posed. It is also clear that this singularity cannot 
have any physical consequences. It arises from flips of triggering spins which are very far 
away from the central spin, and these flips cannot affect the dynamical behavior of the 
central spin since they lead to miniscule changes in the bias seen by the latter. Fortunately, 
it is not necessary to take any special precautions about this singularity in the numerical 
integration, since it is automatically regulated by the finiteness of the sample. 

We can understand the last point further as follows. In the collision terms, the integrals 
over e' are equivalent to a sum over sites. Consider, for example, the two-spin loss term 
(process I). We can write 



dfa 
dt 



= -/*(*) £?W I d -^g{\{e-e>y) 



(4.23) 
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Now, 



/ Y^-^y) = E/ de'8(e'-(e-2K ij y) 



E 1 - ( 4 - 24 ) 



In the same way, for the gain term (process II), we have 



^^(e + 2^). (4.25) 

3& 



Hence the combined contribution of these two processes can be written as 



dfa 
dt 



Let us suppose that we cut off the sum so that site j lies inside a large sphere (centered at 
site i) containing N K sites. For spins outside this sphere, is very small in magnitude, 
and we may approximate 

{Me + 2a'K l3 ) - f a (e)) = ^ 2 d -^a'K lv (4.27) 

3>N K j>N k 

which vanishes if we divide the sum into subsums carried out over sets of sites related by 
cubic symmetry. 



We also note that the form (4.20) for g(K) is a poor estimate for near neighbor or short- 
distance couplings. It is these couplings which are in the end responsible for the peaked 
form of the bias distribution function, and must therefore be treated correctly taking the 
discontinuous delta-function nature into account. In other words, we must handle the part 
of the integral where e' — e is large as a sum over neighboring sites, and a continuum form 
for g(K) cannot be used. 



B. Numerical Integration 

To integrate the kinetic equations, the bias distribution functions /+(e) and /-(e) were 
approximated by histograms with a scaled bias range [— e max , £ ma x], divided into Nf, bins. As 
in the Monte Carlo simulations, £ max was set to 25. At first sight it appears that we should 
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choose the bin width, Wb = 2e max /Nb, to be much less than w. Since realistic values of w are 
quite small, however, this would require the number of bins, Nb to be rather large. Now a 
larger number of bins leads to a smaller number of spins in each bin, and since we want the 
relative change in this number per time step to be small for accurate integration, it requires 
smaller integration timesteps. Hence too small a bin width is numerically expensive. In 
practice we find that changing Wb only affects how very small changes in bias are handled, 
and it is perfectly acceptable to let the bin width equal the reversible region width 2w, leading 
to N b in the range 10 3 -10 4 . Smaller values of Wb did not lead to appreciably different results. 

Since the kernel is a property solely of the lattice, we determine it once and for all 
before doing any integration as soon as we have decided upon the bin width Wb- All values 
of — Kij are determined for a large sphere of Nk spins with the site i at the center. The 
quantity —K^ is half the amount by which the bias will change at the central site if a spin 
(jj flips from up to down. We therefore compile a histogram from this data with a bin width 
matching that of the bias distribution functions, giving us g((e — e')/2jde. It should be 
noted that it does no good to bin g (K) more finely. Let us denote this bin width by w g . 
The value in a given bin of the g{K) histogram essentially gives us the number of sites that 
can shift the bias by 2K. The value of g (K) in the next bin will gives us the sites that shift 
the bias by 2(K + w g ). If 2w g is smaller than Wb, we are needlessly differentiating between 
sites that have the same effect as far the evolution of the histogram for f a {e) is concerned. 
We therefore choose 2w g equal to Wb- 



Combining Eqs. (4.23) and (4.24) for the contribution to df/dt from process I (the loss 
term), we get 

-Nk(T\i + T 1 j)/ cr (e) (4.28) 



dfa 
dt 



We can combine this with the term from process II to write 



dfa 

dt 



= \Ys T -°'°' I de'g R (\{e'-e)a')U{e') 
i+ii z „, J 



(4.29) 



where 

g R (K) = g(K) - N K 6(K) (4.30) 
is a regulated density of couplings. 



The form (4.29) once again shows that the infrared divergence mentioned in the previous 
subsection is a numerical nonissue. Once we have settled on a bin width Wb/2 for the 
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histogram of g{K), there is no point in increasing Nk beyond a certain value. Beyond that 
value, we only change the number in the central bin around K — 0. The delta-function 



term in Eq. (4.30) goes into the same bin, so Qr{K) does not change. Thus the numerics 
are essentially done for an infinite system, and the width of the bin serves as a cutoff that 
effects the principal value integral. 

It also pays to define the population transfer rate 

T(e,e')= 1 -J2 T ^9n(W-^y)- ( 4 - 31 ) 

a 1 

When we multiply this quantity by the bin width u> 6 , the quantity §<7ii(§(e' — e)a')wb is the 
number of sites capable of triggering a shift in bias from e' to e, and T& a i is the spin-flip rate 
averaged over all sites, so WbT(e,e') is the rate at which population is transferred from the 
bin containing the bias e' to the bin containing e. It also gives the entire collision integral a 
nifty matrix multiplication form, 



dfg 
dt 



= [ deT(e,e')/ ff (e') 
i+ii J 



(4.32) 



and the complete kinetic equation can be written as 
df a dfa de a 



T, a (e)f a (e)+T a ,(e)Me) + j de'T(e, e')/ CT (e'). 



dt de dt 

The actual integration of the kinetic equations is simple once the histograms for / + , /_, 
and gn(K) are set up. In each timestep, there is some exchange of the populations in the 
central bins of the / + and /_ histograms, corresponding to the reversible region, as dictated 



by Eq. (4.3). The reversible region is swept from — e max to e max for LZS runs, while it 
is allowed to remain static at the origin for relaxation runs. Next, the collision terms are 



evaluated for every bin in the two histograms, given by Eq. (4.29 ), and the integration moves 
to the next timestep. 



C. Kinetic Equation Results 

As a test of the kinetic equations and our numerical integration procedure, we first applied 
them to the problem of magnetic relaxation with zero external bias. The equations were 
integrated up to a time lOr where r is the characteristic time for relaxation, given by 

r = j!*. (4.33) 
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FIG. 7: (Color online) Comparison of Monte Carlo (MC) and kinetic equation (KE) results for 
magnetic relaxation with w = 0.05. MC data was obtained using ./V = 82,519 and averaged over 
20 runs. KE data was obtained using 2,001 bins and a scaled bias range (—50, 50). 

(Note that this time scale is very long in an absolute sense since it varies as A -2 , and A, 
being a tunnel splitting between deep levels on opposite sides of the anisotropy barrier, is 
very small.) The resulting demagnetization curve is plotted in Fig. [7] along with the results 
of the Monte Carlo simulation. As can be seen, it agrees extremely well with the Monte 
Carlo simulation for short times. However, after t/r > 1, the two curves begin to diverge, 
with the kinetic equations giving a higher value for the magnetization than the Monte Carlo 
simulation. This difference can be understood by comparing the bias distributions given by 
the two methods. These distributions are shown in Fig. |HJ 

For short times, t < O.lr, the bias distributions match fairly well, with each peak in the 
Monte Carlo data also present in the solution to the kinetic equations. However, after t ~ r, 
the bias distribution from the Monte Carlo is bounded by about ±10E dm , and still shows 
sharp side peaks, while that from the kinetic equations is significantly more spread out and 
lacking the peaks. This difference is due to the fact that the kinetic equations do not carry 
information about the orientations and biases of a spin's neighbors while the Monte Carlo 
simulation does. More specifically, there is a short distance correlation between biases, as 
we now discuss. Consider a central spin and its 6 nearest neighbors all initially in the up 
state and with zero bias. If one of the neighboring spins flips, the bias on the central spin 
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FIG. 8: Comparisons of how the bias distributions evolve during magnetic relaxation as found 
by Monte Carlo (MC) simulations and by solving the kinetic equation (KE). Left and right hand 
panels: KE and MC for t/r = 0.1 (top) and 1.0 (bottom). The scaled reversibility region width 
is w = .05. MC data was obtained using N = 82,519 and averaged over 20 runs. KE data was 
obtained using 2,001 bins and a scaled bias range (—50, 50). 



will change by an amount 8Ed m or —AEdm, moving it out of the reversible region to one of 
the peaks that can be seen in Fig. [8b. Similarly, the bias at the 5 other neighboring sites 
will also change and these spins will also be removed from the reversible region. As the 
sample continues to relax, the bias for the central spin will not have a large probability to 
change significantly because of the fact that 5 of its nearest neighbors are displaced far from 
the reversible region. The sixth spin, the one that originally nipped, is still sitting in its 
original near-zero bias, and therefore has a higher probability to flip back, but this would 
cause the bias for the central spin to move back closer to the origin. Thus, by considering a 
spin's neighbors, we can see that once a spin acquires a large bias, there is a low probability 
for the bias on that spin to change yet again by a large amount, and if such a change does 
occur it is more likely to make the bias regress back toward the mean. This explains why the 
bias distribution in the Monte Carlo simulation is more bounded and peaked. The kinetic 
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equations do not keep track of these correlations, and thus yield too broad a distribution. 

We then applied the kinetic equations to the problem of a swept field. In this case, the 
equations are much more successful, and, as shown in Fig. |9j we find that the results are 
in good agreement with the Monte Carlo simulation for values of the scaled sweep rate, v, 
as low as 5, where the Kayanuma model does quite poorly. However, the limitations of the 
kinetic equations again show up for v < 5. 




FIG. 9: (Color online) Mf vs. v as given by the Kayanuma model, our Monte Carlo data, and the 
kinetic equation. For the latter two we used a width parameter of w = 0.05 and for the kinetic 
equations, we used a bin width of Ed m /S0. 



V. CONCLUSIONS 



We have carried out Monte Carlo simulations of swept-field experiments on molecular 
magnetic solids based on the microscopic picture of spin reversal developed in Refs. [T5l 
[3i] . We find that these simulations provide a very good picture of the time evolution 
of the entire system, and agree fairy well with experiments quantitatively. In order to 
understand the simulations, we have also developed a kinetic equation for the distribution 
of single-site spin and bias distribution. This kinetic equation also provides a quantitatively 
accurate description of experimental data even for quite low sweep velocities. However, the 
kinetic equation fails at very low sweep velocities, since it is then incapable of accounting 
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for important short-distance bias correlations. Expanding the kinetic equation approach to 
include two-site distributions is non trivial and difficult. Nevertheless, the kinetic equation 
should be capable of describing many more experiment protocols, and we hope to do this in 
the future. 
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